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Abstract 

The aim of this paper is to compare three regularized particle filters in an online data 
processing context. We carry out the comparison in terms of hidden states filtering and 
parameters estimation, considering a Bayesian paradigm and a univariate stochastic volatil- 
ity model. We discuss the use of an improper prior distribution in the initialization of the 
filtering procedure and show that the regularized Auxiliary Particle Filter (APF) outper- 
forms the regularized Sequential Importance Sampling (SIS) and the regularized Sampling 
Importance Resampling (SIR). 

Keywords: Online data processing; Bayesian estimation; regularized particle filters; stochas- 
tic volatility model 
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1 Introduction 

The analysis of phenomena, which evolve over time is a common problem to many fields like 
engineering, physics, biology, statistics, economics and finance. A time varying system can be 
represented through a dynamic model, which is constituted by an observable component and an 
unobservable internal state. The hidden states (or latent variables) represent the informations 
we want to extrapolate from the observations. 

In time series analysis, ma ny approaches have been used for the e stima tion of dynamics 
models. The seminal works of Kalmanl ( 1960l ) and Kalman and Bucv ( 1960l ) introduce filter- 
ing tech niques (the Kalm an-Bucy filter) for continuous valued, linear and Gaussian dynamic 
systems. iMavbeckl (|1982i ) motivates the use of stochastic dynamic systems in engineering and 
examines the estimation probl ems for s tate s pace models, in both a continuous and a discrete 
time framework. In economics, iHarvev! (1 19891 ) studies the state space representation o f dynamic 
structural models and uses Kalman filter for hidden states filtering. Hamilton] ( 19891 ) analyzes 
nonlinear time series models and introduces a filter (Hamilton-Kitagawa filter) for discrete time 
and discrete valued dynamic systems with a finite number of states. 

In this pap er, the online data p rocessing problem is considered. In these situations, as 
pointed out by iLiu and Chenl ( 19981 ) . Markov Chain Monte Carlo (MCMC) samplers are too 
much time demanding. To overcome this difficulty, some sequential Monte Carlo techniques have 
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been recently developed. iDoucet et al,l (|200ll ) provide the state of the art on these methods. 
They discuss both applications and theoretical convergence of the algorithms. 

The contribution of this work is the comparison of three types of regularized particle fil- 
ters - the regularized Sequential Importance Sampling (SIS), the regularized Sampling Im- 
portance Resampling (SIR) and the regularized Auxiliary Particle Filter (APF) - when the 
model parameter s are unknown. T he online estimation of m o del paramet e rs is a difficult task 
(iKitagawal (Il998h: IStorvik! (l2002ll: berzuini and Gilkd (EoOlll: iFearnheadl (l2002h: iDiuric et al ' 



2002yStorvikl(j20n2h : lAndrieu and Doucetl (j2003h : lDoucet and Tadic 



We consider here the Bayes i an paradigm and the regularization (see 



(2003); 



Poison et al 



Chen and Havkinl (|2002j)) 



(2002)) 



Rossil (|2004l ) based on a 



approach of lOudianel (j2000h : iLiu and WestJ toll ): IMusso et al.1 (j200l1 i 
kernel approximation in the parameter-augmented state space. We also discuss the initialization 
of the filtering procedure. 

This work is structured as follow. Section [2] introduces the general representation of a 
Bayesian dynamic model and presents the stochastic volatility model. Section reviews some 
regularized particle filters. Finally, Section H] presents the results. 



2 Bayesian Dynamic Models 



We introduce the general formulation of a Bayesian dynamic model and show some fundamental 
relations for Bayesian inference on it. O u r defi ni tion of dy n amic m odel is general en o ugh t o 
include the models analyzed in iKalmanl (Il960h. iHamiltonl (|l994l ). ICarter and Kohnl (jl994T ). 
Harrison and West ( 19891 ) and in Doucet et al.l ( 2001 ). Throughout this work, we use a n otation 
similar to that one commonly used in particle S^s literature (see boucet et all ». 

We denote by {x t ; t G N}, x t G X C R n * , the hidden states of the system, by {y t ; t G No}, 
yt G y C M. n v, the observable variables and by {6 t ; t G N}, ( £ 6 C R ne , the parameters of the 
model. We denote by xo : j = (xo, . . . ,x t ) the collection of hidden states up to time t and with 
x_j = (xo, . . . ,xt_i,xt + i, . . . ,xr) the collection of all hidden states without the t-th element. 
We use the same notations for the observable variables and parameters. 

The Bayesian state space representation of a dynamic model is given by: 



yt ~ p(yt|x t ,0 t ,yi :t _i) 

(x 4 ,0 t ) ~ p(X-t,Ot\xO:t-l,Oo:t-l,yi:t-l] 

xo ~ p(xol^o) 

O ~ 7r(0 O ) 



measurement density , 
transition density , 
initial density , 
prior density , 



for t = 1,...,T. 

We suppose that p(x t , 0t|x O: £_i, O :t-i, yi:t-i) = p(*t, Ot\x-t-\, Ot-u Yi-.t-i)- We also assume 
that the parameters are constant over time: the transition density of the parameters is then 
5ff t _ 1 (6 t ) with initial value Oq = 9, 6 x (y) denotes the Dirac's mass centered in x. 
In that case, the joint transition of hidden states and parameters is: 

p(x t , 0t|x*-i, Ot-i, yv.t-i) = p(xt|x 4 _i, e t , yi-.t-i)Se t ^ (Qt) ■ 

Let us denote by zj = (xj,0t) the parameter-augmented state vector and by Z the corre- 
sponding augmented state space. For such models, we are interested in the prediction, filtering 
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and smoothing densities which are given by: 



P(z*+l|yi:t) = J P(xt+l\xt,0t+l,yi:t)50 t (0t+l)p(zt\yi:t)dz t , (1) 

p(yt+i|yi:t) = J p(yt+i\zt+i,yi:t)p(zt+i\yi:t)dz t+ i , 

I i \ P(yi+l|Zi+l,yi:i)p(Zi+l|yi:i) fo s 
P{Zt+l\yi:t+l) = j i r , \4) 

p(y*+i|yi:t) 

P{Zs\yi:t) = P{Z s \yi:s) / j v dZ a+1 , 3<t. 

Jz p(z s +i|yi :s ) 

Due to the high number of integrals that must be solved, previous densities may be difficult to 
evaluate with general dynamics. Some Monte Carlo simulation methods, such as particle filters, 
allow us to overcome these difficulties. 

As an example, let us consider the stochastic volatility model. Two of the main features 
of the financial time series are time varying volatility and clustering phenomena in volatility. 
Stochastic volatility models widely used in finance have been introduced, in order to account 
for these features. Let yt be the observable variable with time varying volatility and xt the 
stochastic log-volatility process. An example of stochastic volatility model is: 

y t \x t ~ Af(0,e xt ) 
x t \x t -i,0 ~ N (a + 4>x t -i,a 2 ) 

x \e ~ aa(o, C t 2 /(i-0 2 )) 
e ~ tt(0) 

where = (a, log((l + </>)/(l — (/>)), log(cr 2 )). The choice of 7r(0) will be discussed in Section HI 
Fig. [T] shows two simulated paths of yt and xt- 

In the next section, we deal with the problem of parameters and states joint estimation in 
a kernel-regularized sequential Monte Carlo framework. 



3 Regularized particular filters 



For making inference on the Bayesian dynamic model given in Section [5] in an online data 
processing context, MCMC algorithms are too much time demanding. Sequential impor- 
tance sampling and more advanced sequential Monte Carlo algorithms called Particle Filters 
(jDoucet et all l200ll ) represent a promising alternative. The main advantage in using particle 
filters is that they can deal with nonlinear models and non-Gaussian innovations. In contrast 
to Hidden Markov Model filters, which work on a state space discretized to a fixed grid, par- 
ticle filters focus sequentially on the higher density regions of the state space. This feature is 
comm on to one of the ea rly sequential methods, the Adaptive Importance Sampling algorithm 
due to lWestl t|l992l . Il993h . 

Different particle filters exist in the literature and different simulation approaches like rejec- 
tion sampling, MCMC and importance sampling, can be used for the construction of a particle 
filter. In this work, we present some kernel-regularized particle filters, which combine the im- 
portance sampling reasoning with a suitable modificatio n of the importance weig hts. The reg- 
ularization approach we use is the same than the one of iLiu and West (|200lh andlMusso et al 



(|200lh . This approach relies upon a kernel-based reconstruction of the empirical filtering den- 
sities which produces a systematic modification of the true importance weights. 



3 



Daily 




to 



200 400 600 800 1000 

Time 

Weekly 




200 400 600 800 1000 

Time 



Figure 1: Simulation paths for x t {grey line) and yt {black line). Upper plot: daily dataset 
[a = 0, cf) = 0.99 and a 2 = 0.01). Bottom plot: weekly dataset {a = 0, <j) = 0.9 and a 2 = 0.1). 



3.1 Regularized SIS 

Let us start from the non-regularized SIS. We assume that at iteration t > a properly weighted 
particle set {x£, dttli}iLii approximating the filtering density p(x t , d t \yi-t), is available. The 
empirical distribution corresponding to this approximation is: 

N 

PN(xt,0t\yi*) = ^2il\ x i,9i)( x t,o t ) . (3) 

i=l 

The particles set, {x£, l t , ~i\}f = i, can be viewed as a random discretisation of the state space 
X x O with associated probability weights {jl}iLi- Thanks to this discretisation, it is possible 
to approximate the prediction and filtering densities given in (JTJ) and @: 

N 



i=i 

N 

pjv(x t+ i,0 m |y 1:t+ i) oc ^7^>(y t+ i|x t+ i,0 t+1 ,yi : t)p(xi + i|xJ,0 t +i,yi:i)^}(0t+i> 



i=i 



The goal is now to obtain iV particles {xj +1 , 0j +1 , 7^ , from the filtering density in (|2|). 
It is proposed to sample (xj +1 , 6\ +l ) according to the importance density g(- |x|, 0\, yi : t+i). The 
importance weight of particle (x£ +1 ,0J +1 ) is then calculated using the recursive formula: 

t+1 * 9(xj +1 ,6»j+i|xj,6>j,y 1:m ) 
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The choice of an optimal importance density g(-|x|, 0\, yi;t4-i), that is, a density which min- 
imizes the variance of th e importance weights is discussed in IPitt and ShephardI dl99^ and 
Crisan and Doucet ( 200d ). In many cases, it is not possible to use this optimal importance 
density as the weight updating associated to the this density does not admit a closed-form 
expression. In that case, the transition density of the parameter-augmented state vector repre- 
sents a natural alternative for the importance density. Indeed, the transition density represents 
a sort of prior at time t for the parameter-augmented state vector (xj +1 , 0\ +1 ). 

In our case, due to the presence of the Dirac point mass at the numerator of the weights it 
is impossible to modify over the filtering iterations the particle values for the parameters. In 
practice due to the loss of particle diversity in the parameter space, the weights will tend to 
zeros and of course stay zero for ever, so we are facing a problem of degeneracy of the empirical 
filtering distribution. This scenario motivates particle filtering methods known as regularized 
particle filters. In order to avoid the degeneracy problem and to force the expl oration of the 
param eter space toward region s which are not covered by the prior distribution, iLiu and West 
(|200lh and lMusso etail (|2f)0lh ; propose to use a regularized version of the filtering density. This 
approach results in the modification of the weights in @ and the definition of a new set of 
weights: 



t+n x tj^t+i,yi:t)^/i(^+i -Q\ 



q(-xi +1 ,e l t+1 \4,ei, yi:t+1 ) 

where Kh(y) = h~ d K(y /h) is a regularization kernel, K being a positive function defined on 
W 16 and h a positive smoothing factor (bandwidth). 

The modification of the importance weights defined in (|3]) results from two steps. The first 
one is the regularization of the empirical density in @ by a kernel estimator: 



N 



PN(*t,Ot\yi:t) 



i=i 



The second one is the application of an importance sampling argument to the approximated 
filtering density: 



N 



p$(x m ,0 m |yi :t+ i) = ^2u l t p(y t+1 \x t+1 , m ,yi :t )p(x m |xj, O t+1 ,y 1:t )K h (O t+1 - 0\) 



i=l 



Th e con vergences result s associated with this type of approximation are recalled in lMusso et al 
toil ) and lOudiand (2000). Under some usual conditions on the kernel, when the number of 
particles increases to infinity, the regularized empirical density converges to the right one for 



various criteria. For instance, we have pjy — >£i p. 

Thanks to this approximation, the regularization kernel becomes the natural choice for the 
parameters proposal distribution. Thus, we sample (x£ +1 ,#J +1 ) according to: 

q(xt+i\*i,0t+i,yi:t+i)K h (0t+i ~ 0\) . 



In that case, we have: 



i p(ym 1*4+1 , o\+i » yi:tMxt+i |*t , 0t+i , yi:t) 



\*l, e t+i,yi:t+i) 

In Algorithm [H we give a pseudo-code representation of this method. 
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Algorithm 1 - Regularized SIS Particle Filter - 

■ At time t = 0, for i = 1,...,N, simulate Zq ~ p(zo) an d set ujq = 1/N 

■ At time t > 0, given {x\,6\,uj l t }f =l , for i = 1, ... ,N: 

1. Simulate 0\ +1 ~ K h (0 t+1 - 0\) 

2. Simulate xj +1 ~ <?(x t+1 |x|, l t+1 , yi :t +i) 

p(yt+i l x m , e \+i , yi:t)p(*l + i |x|, , yi :t ) 



5". Update the weights: oc 



g(xj +1 |xj,^ +1 ,yi : i + i; 



3.2 Regularized SIR 

As it is well known in the literature (see for example Arulampalam et al. ( 200ll )). basic SIS 



algorithms have a degeneracy problem. After some iterations the empirical distribution degen- 
erates into a Dirac's mass on a single particle. T his due to the fact t hat the variance of the 
importance weights i s non-decreasing ove r time (see boucet et ai](|2000l )). In order to solve this 



degeneracy problem, Gordon et al. ( 19931 ) introduce the SIR algorithm. This algorithm belongs 



to a wider class of bootstrap filters. At each iteration, a resampling step is used to generate 
a new set of particles. After this resampling step, the weights of the resampled particles are 
uniformly distributed over the particle set. 

In the initial SIR, the resampling step is done at each iteration o f the algorithm . Th is 
systematic resampling can introduce extra Monte Carlo variations, see iLiu and Chenl dl998l ). 



This can be reduced be doing resampling only when the Effective Sample Size (ESS) is small. 
The ESS measures the overall efficiency of an importance sampling algorithm. The ESS is a 
function of the coefficient of variation of the importance weights. At iteration t, the empirical 
EES is 

ESS t = — 



N 



1 + 

1=1 

In Algorithm [21 we give a pseudo-code representation of this method. 




Algorithm 2 - Regularized SIR Particle Filter - 

■ At time t = 0, for i = 1, . . . , N , simulate Zq ~ p(zo) an d set uJq = 1/N 

■ At time t > 0, given {x\, 6\, u)\}f =1 , for i = 1, . . . , N : 

1. Simulate 0\ +1 ~ K h (G t +i - 0\) 

2. Simulate xj +1 ~ g(x t+ i|x|, d\ +1 , yy. t +x) 



3. Update the weights: <-o^+i oc uj, 



f p(y t+ ilxj +1 ,g f+1 ,yi : t)p(x| +1 |x|,fl m ,yi :t ) 
g(x| +1 |x|,^ +1 ,y 1:m ) 



I IfESSt+i < k, simulate {x^,^}^ from H+iA + lA+i\i=i (Multi- 
nomial resampling) and set = 1/N. 



The value of k < N is calibrated depending on the problem. 
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3.3 Regularized APF 

Due to the resampling step, the basic SIR algorithm produces a progressive impoverishment (loss 
of diversity) of the information contained in the particle set. To overcome t his difficulty, many 
soluti ons have been proposed in the literature. We re fer to the APF due tolPitt and Shephard 
and to the regularized APF algorithm due to iLiu and Westl (120011 ) . In order to avoid 
the resampling step, the APFs use the particle index (auxiliary variable) to select most rep- 
resentative particles in the proposal of the new particles. The regularized joint distribution of 
parameter-augmented state vector and the particle index is: 

p$(xt+i, e t+ x, *|yi:i+i) oc p(y t+1 |x t+1 , 0£ + i,yi :t )p(x t+ i|x£, l t ,y 1: t)K h (0t+i - 0\)uj l t . 

A sample approximating that distribution can be obtained by using the proposal: 

q(4+l , Ol+l , f |yi:t+l) = P(*l+l \*t\ , yv.t)Kh{0\ +l - 0{)q{f \yi:t+i) 

where 



g(/|yi:t+i ) « p(yt+i\(4+i > m m > yi^ ) w \ 



t i 



At A1 

H t+ i and m t , 1 are evaluated using the initial particle set. Therefore, the importance weight of 
particle (xj +1 , 0\ +1 , f) is: 

p(y t +i\xl + i,Ol +1 ,yi:t) 



p(yt+i|^ +1 ,m^ +1 ,yi :t ) 
In Algorithm [3] we give a pseudo-code representation of the regularized APF. 



Algorithm 3 - Regularized Auxiliary Particle Filter - 

■ At time t = 0, for i = 1, . . . ,N, simulate z l ~ p(zo) an d set ujq = 1/N 

■ At time t > 0, given {xj, 6\,oj % t }f =1 , for i = 1, . . . ,N: 

1. Simulate f ~ l(j\yv.t+l) with j € {1, . . . , N} (Multinomial sampling) 
where /i J t+1 = E(xt+i|x^, 0f) and m? t+l = E(0( + i|0^) 

2. Simulate 0\ +l ~ K h (6 t+1 - 0{) 

3. Simulate x£ +1 ~p(x t+ i|x| ,^ +1 ,yi : t) 

p(yt+i|x|+i»^t+i»yi:t) 



4- Update particles weights: oj\ + i oc 

p(yt+i\rft +1 ,ml' +1 ,yi.. t ) 



We can say that, in the APF, the selection step is done before simulating the hidden states. 
This selection depends on the current value of the observable. Therefore, 

• the APF is a standard way to construct a proposal distribution for the hidden states that 
depends on the current of the particle; 

• as we will see after, to use this selection step and the transition distribution as proposal 
distribution for the hidden states, results in a good proposal distribution. 

In the next section, we compare the performances of some regularized SIS, SIR and APF 
for the stochastic volatility model. 
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4 Application to the stochastic volatility model and comparison 



In this section, we apply the three regularized particle filters to the stochastic volatility model 
presented in Section [21 We assume that the initial value of the SV process follows the stationary 
distribution: 

x ~Af(0,a 2 /(l-4> 2 )). 
For the parameters /3, a and <j), we assume the prior 

p(/3 2 ,<A, f T 2 ) = l/(a/?)I ( _ ljl) (0), 

where (5 = e a . We constrain the parameter </> to take values in the open interval (—1,1) in 
order to impose the usual stationarity condition. As we use an improper prior, it is not possible 
to use the prior distribution for initializing the three particle filters. We need to start with a 
proper weighted sample {x l , Gq, WoH=i- We propose to: 

1) start the sequential filtering procedure at least at the value n of t, such that the posterior 
distribution of the parameters given all the observations up to time n is well defined: for 
the considered SV model, this corresponds to set n > 2; 

2) use a Markov Chain Monte Carlo (MCMC) algorithm to create a sample with uniform 
weights. 

For the prior given above, the full conditional distributions are 
(3 2 \--- ~ jg^y 2 exp(-x t )/2,(n-l)/2^ , 

a 2 \--- ~ 1G (^(x t -0s^i) 2 /2 + x?(l-0 2 ),(n-l)/2) , 



a=2 

/ n-l 



tt(4>\---) oc (l-0 2 ) 1/2 exp l-(f> 2 ^2x 2 -2(f>^2xtxt-i \ /2a 2 I 



(-i,+i) \<P) i 



t=2 t=2 



n(x t \ ■■■) oc exp _ a ~ fat-i) 2 - (xt+i - a- 4>x t ) 2 ) - ^ (x t + y 2 exp(-x f )) 

The full conditional distributions of and %% are not conventional and the standard Gibbs 
do es not apply. We p ropose to use the Metropolis-Hastings within Gibbs algorithm studied 
Celeux et al. ( 20061) . A detail ed description of the proposal distributions for <fi and xt can 



m i^eieux et ai.l (izuua). A aetaiicci description 01 tne proposal aistriDutions lor <p ana xt can 
be found in ICeleux et all l|200d ). In that paper, the authors compare this MCMC scheme 



to an iterated importance sampling one. Note that one could alternatively use this iterated 
importance sampling algorithm to create a first weighted sample. 

Given the initial weighted random sample {x l t , 6\, w|}. =1 , where t = (at, log((l + <^t)/(l — 
<f>t)), log(o" 2 )), if we use the transition density as proposal distribution for the hidden states, the 
regularized SIS performs the following steps: 



8 



SIS 




SIS 




200 400 600 800 
SIR-r-p 




200 400 600 800 
SIR-r-p 




200 400 600 800 
APF 




200 400 600 800 
APF 




200 400 600 800 



200 400 600 800 



Figure 2: Daily (left column) and weekly (right column) true (black line) and filtered (grey line) 
log-volatility. 



For n < t < T — 1 and for i = 1, . . . , N : 

(i) Simulate 0\ +1 ~ M (a0\ + (1 - a)6 t , b 2 V t ) where V t and t are the 

empirical covariance matrix and the empirical mean respectively, 

a G [0, 1] and b 2 = (1 - a 2 ) , 



(ii) Simulate x\ +l ~ M (aj +1 + 4>\ +1 x\, (o- 2 )* +1 J , 
(Hi) Update the weights as follow 



'A ex P | -\ [yf+i ex P {- x l+i) + x l+i) } • 



In the following, we call the previous scheme SIS. Instead of the transition density, we can use 
a proposal distribution which depends on the current value of the observable. For instance, 
we can resort to the proposal distribution used fo r the hidden states, in the MCMC initializa- 
tion step. This proposal has been introduced by IShephard and Pittl (119971) . It is based on a 
qu adratic Taylo r expa nsion of exp(xj), more details can be found in IShephard and Pittl (|l997h 
or Celeux et al. It defines a new SIS which is called SIS-p in the following. 

The results of a typical run of the regularized SIS on the synthetic dataset in Fig. [H with 
N = 10, 000 particles and n = 100 for the Gibbs initialization, are given in Fig. from[2]to[6l We 
can see (last row in Fig. [2]) that after a few iterations the filtered log-volatility does not fit well 
to the true log-volatility. We measure sequentially the filtering performance of the regularized 
SIS by evaluating the cumulated Root Mean Square Error (RMSE). It measures the distance 
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between the true and the filtered states and is defined as: RMSE t = {jYlu^i^u — z u ) 2 }2, 
where &t is the filtered state, which includes also the parameter sequential estimate. The RMSEs 
cumulate rapidly over time in both daily and weekly datasets (see upper and bottom plots in 
Fig. [3|). The poor performance of the regularized SIS is due to the fact that the empirical 
posterior of the states and parameters degenerates into a Dirac's mass after a few iterations. 
The ESSs in Fig. 0] show that the regularized SIS degenerates after 30 iterations in both the 
daily and weekly datasets. We give some results on SIS-p in the following. 

If we use the transition density as proposal distribution for the hidden states, the regularized 
SIR performs the following step: 

For n < t < T — 1 and for i = 1, . . . , N: 

(i) Simulate 0\ +1 ~N (a0\ + (1 - a)6 t , b 2 V t ) where V t and 6 t are the 

empirical covariance matrix and the empirical mean respectively 
and a G [0, 1] and b 2 = (1 - a 2 ) , 

(ii) Simulate x\ +1 ~ M (a\ +1 + 4>\ +l x\, {<? 2 )\ + -^) > 
(Hi) Update the weights 

w\ +l oc w\ exp | -i [y 2 +l exp{-x l t+l } + x l t+1 ] j , 

(v) If ESS m < k, simulate z\ +1 ~ 12j=i w t+l^ x p K+i) and set w t+i = 
l/N. 



If k = N, the resampling step is done all the time. In that case, we call SIR the previous scheme. 
After some numerical experiments, we have found that a good value for k is k = 0.9 x N. In 
that case, the resampling step is done at regular time intervals and w e called SIR-r the re s ulting 



algorithm. Moreover, as for the SIS, we can resort to the proposal of lShephard and Pittl (| 199 71 ) 
for the hidden states. In that case, with n = N, the corresponding algorithm called SIR-p. 
With k = 0.9 x N the corresponding algorithm is called SIR-r-p. 
The regularized APF performs the following steps: 

For n < t < T — 1 and for i = 1, . . . , N: 

(i) Simulate f ~ q(j) cx £f=i w$Af(y t+ i\4 +1 )6 k (j) where = + a\ , 



(ii) Simulate 0\ +l ~N \a6{ + (1 - a)O t ,b 2 V t j where V t and 6 t are the 
empirical variance matrix and the empirical mean respectively 

and a G [0, 1] and b 2 = (1 - a 2 ) , 



(Hi) Simulate x\ +l ~ M [x t+1 \a\ +l + 4>\+\A > {v 2 )]^) > 
(iv) Update the weights 



w\ +1 oc exp |~ y 2 t+l (exp{-x l t+l } - expj-^j) + x\ +l - \ 
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Figure 3: Daily (upper plot) and weekly (bottom plot) Root Mean Square Errors for the regu- 
larized APF (solid line), SIR-r-p (dashed line) and SIS (dotted line) over iterations. 




Figure 4: Daily (left column) and weekly (right column) Effective Sample Sizes over iterations. 
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Figure 5: Evolution on daily dataset of the empirical posterior distributions of a, 4> an d o -2 
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Figure 6: Evolution on weekly dataset of the empirical posterior distributions of a, <j> and a 1 
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Daily Data 


6 


SIS 


SIS-p 


SIR SIR-p 


SIR-r 


SIR-r-p 


APF 


a 


u.uu i 19 


0.00945 


0.00885 0.00925 


0.00315 


0.00912 


0.00065 


4> 


0.66767 


0.83264 


0.12433 0.13355 


0.15252 


0.13456 


0.00855 


a 2 


0.89327 


0.87910 


0.00676 0.00670 


0.00643 


0.00654 


0.00506 



Table 1: Mean Square Errors of the estimators of a, (ft and a 2 . The Mean Square Errors are 
estimated using the last iteration of the 10 independent runs of the filters. 



Weekly Data 


e 


SIS 


SIS-p 


SIR SIR-p 


SIR-r 


SIR-r-p 


APF 


a 


0.00534 


0.00487 


0.00589 0.00442 


0.00380 


0.00431 


0.00016 


<t> 


0.51290 


0.55648 


0.05292 0.03754 


0.04885 


0.03687 


0.00029 


a 2 


0.70540 


0.71242 


0.00010 0.00009 


0.00009 


0.00009 


0.00008 



Table 2: Mean Square Errors of the estimators of a, <j) and a . The Mean Square Errors are 
estimated using the last iteration of the 10 independent runs of the filters. 



Note that, following [Pitt and Shephardl (119991 ) . one could alternatively use in the selection step 
a value of based on the Taylor expansion of the likelihood at time t + 1. For the three 
regularized particle filters, we have used a Gaussian kernel where the parameter a is fixed 
following the usual optimal criterion. 

We apply the regularized SIR-r-p and APF with N = 10, 000 and n = 100 to the weekly 
and daily datasets of Fig. [T] and obtain the results given in Fig. from [2] to [6j The regularized 
SIR-r-p and APF outperform the regularized SIS in terms of ESSs and cumulated RMSEs. 
The ESSs can detect the degeneracy in the particle weights, but is not useful to determine the 
presence of another form of degeneracy, that is the absence of diversity in the particle values. 
The histogram of the empirical filtering distribution allows us to detect this second form of 
degeneracy. 

As our work deals with the sequential estimation of the parameters, we choose to show 
the histogram of the parameters posterior. Fig. [5] and [6] exhibit the evolution over the filters 
iterations of the posterior of the parameters a, 4> and a 2 . In both the daily and the weekly 
datasets, after a few iterations the empirical posterior of the regularized SIR-r-p degenerates 
into a Dirac's mass. 

To confirm the previous results, we have done ten independent runs of the seven algorithms: 
SIS, SIS-p, SIR, SIR-p, SIR-r, SIR-r-p and APF. The weekly and daily datasets vary across 
the 10 experiments. Our simulation study confirms the results of the single-run experiment. 
Fig. [7] and [8] show a comparison between the regularized schemes in terms of RMSEs. The 
RMSEs are estimated over the 10 independent runs of the algorithms. The regularized SIR-r-p 
and APF outperform the others algorithms in both the daily and the weekly datasets. The 
estimated Mean Square Errors for the parameters a, (j) and a 2 (see Table [JJ and [2]) , based on 
10 independent runs of the filters, show that the regularized APF outperforms all the others 
schemes in term of parameters estimation. 

Fig. [9] and [10] show a comparison of the filters in terms of ESS. As one could expect, in all the 
independent runs the regularized SIS and SIS-p weights degenerate after a few iterations. We 
also observe that the ESSs of the regularized SIR, SIR-r and APF are substantially equivalent. 
On the other hand, the average ESSs of the SIR-p and SIR-r-p are lower than the ones of the 
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SIR, SIR-r and APF. 



5 Conclusion 

In this work we bring into action the kernel regularization technique for particle filters and deal 
with the online parameter estimation problem. While the regularized APF has been already 
used for the parameter estimation, the regularized versions of SIS and SIR have not been 
considered to that aim. We focus on the joint estimation of the states and parameters and 
compare some algorithms on a Bayesian nonlinear model: the Bayesian stochastic volatility 
model. As we expected, we find evidence of the degeneracy of two different regularized SIS. 
Finally, we find that, in terms of parameters estimation, the regularized APF outperforms all 
the others schemes. 
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Figure 7: Comparison on daily datasets of average cumulative Root Mean Square Errors between 
the true and the filtered log- volatility (black line). We represent the area between maximum 
and minimum cumulative Root Mean Square Errors (grey area). 
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Figure 8: Comparison on weekly datasets of average cumulative RMSEs between the true and 
the filtered log- volatility (black line). We represent the area between maximum and minimum 
cumulative Root Mean Square Errors {grey area). 
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Figure 9: Comparison on daily datasets of average Effective Sample Sizes (black line). We 
represent the area between maximum and minimum Effective Sample Sizes (grey area). 
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Figure 10: Comparison on weekly datasets of average Effective Sample Sizes (black line). We 
represent the area between maximum and minimum Effective Sample Size (grey area). 
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